home *** CD-ROM | disk | FTP | other *** search
/ United Public Domain Gold 4 / United Public Domain Gold 4.iso / fredfish / ff.0164.dms / ff.0164.adf / Newton / complex.c < prev    next >
C/C++ Source or Header  |  1988-11-22  |  3KB  |  112 lines
  1. /* complex.c:    Complex arithmetic routines.
  2.  *
  3.  * Written by Daniel Barrett.  100% PUBLIC DOMAIN. */
  4.  
  5. #include "decl.h"
  6.  
  7. /* Complex arithmetic functions Add(), Subtract(), Multiply() and Divide()
  8.  * perform as their names indicate.  They perform their operation on their
  9.  * first 2 arguments, and return the result in the third argument. */
  10.  
  11. Add(a, b, sum)
  12. complex a, b, *sum;
  13. {
  14.     sum->n[REAL] = a.n[REAL] + b.n[REAL];
  15.     sum->n[IMAG] = a.n[IMAG] + b.n[IMAG];
  16. }
  17.  
  18.     
  19. Subtract(a, b, diff)
  20. complex a, b, *diff;
  21. {
  22.     diff->n[REAL] = a.n[REAL] - b.n[REAL];
  23.     diff->n[IMAG] = a.n[IMAG] - b.n[IMAG];
  24. }
  25.  
  26.     
  27. Multiply(a, b, prod)
  28. complex a, b, *prod;
  29. {
  30.     prod->n[REAL] = (a.n[REAL] * b.n[REAL]) - (a.n[IMAG] * b.n[IMAG]);
  31.     prod->n[IMAG] = (a.n[REAL] * b.n[IMAG]) + (a.n[IMAG] * b.n[REAL]);
  32. }
  33.  
  34.     
  35. Divide(a, b, quot)
  36. complex a, b, *quot;
  37. {
  38.     double denom;
  39.  
  40.     denom = SQR(b.n[REAL]) + SQR(b.n[IMAG]);
  41.     if (denom == 0.0)
  42.         printf("Divide by zero error!\n"), exit(20);
  43.  
  44.     quot->n[REAL] = ((a.n[REAL] * b.n[REAL]) + (a.n[IMAG] * b.n[IMAG]))
  45.             / denom;
  46.     quot->n[IMAG] = ((a.n[IMAG] * b.n[REAL]) - (a.n[REAL] * b.n[IMAG]))
  47.             / denom;
  48. }
  49.  
  50.     
  51. SynthDivide(poly, comp, stop)
  52. /* Computes the synthetic division of the polynomial poly[] by (z-comp), 
  53.  * where "z" is assumed the complex variable in our polynomial. */
  54. complex poly[], comp;
  55. int stop;
  56. {
  57.     int i;
  58.     complex prod;
  59.  
  60.     for (i=1; i<=stop; i++) {
  61.         Multiply(poly[i-1], comp, &prod);
  62.         Add(poly[i], prod, &poly[i]);
  63.     }
  64. }
  65.  
  66.     
  67. Iterate(poly, zOld, n, zNew)
  68. /* One iteration of Newton's method.  "zOld" is the old guess of the value
  69.  * of a root of poly[].  "zNew" is the newly calculated guess. */
  70. complex poly[], zOld;
  71. int n;
  72. complex *zNew;
  73. {
  74.     complex quotient;
  75.     Divide(poly[n], poly[n-1], "ient);
  76.     Subtract(zOld, quotient, zNew);
  77. }
  78.  
  79.     
  80. Guess(z)
  81. /* A random complex number, 50% chance of being negative. */
  82. complex *z;
  83. {
  84. #ifdef UNIX
  85. #define    ran    drand48
  86. #endif
  87.     double ran();
  88.  
  89.     z->n[REAL] = ran();
  90.     z->n[IMAG] = ran();
  91.     z->n[REAL] = (ran() < 0.50) ? z->n[REAL] : -(z->n[REAL]);
  92.     z->n[IMAG] = (ran() < 0.50) ? z->n[IMAG] : -(z->n[IMAG]);
  93. }
  94.  
  95.     
  96. BOOL ErrorTooBig(y, z, eps)
  97. /* Compute the Euclidean distance between y and z in the complex plane.
  98.  * This is just our good old friend, the "distance formula".  (Add the
  99.  * squares of y and z, and take the square root.)  Is the distance less
  100.  * than epsilon?  If so, the error is allowable; else, it's too big. */
  101. complex y, z;
  102. double eps;
  103. {
  104.     complex diff;
  105.     double sqrt();
  106.  
  107.     Subtract(y, z, &diff);
  108.     return(    sqrt(SQR(diff.n[REAL]) + SQR(diff.n[IMAG])) < eps
  109.         ? FALSE
  110.         : TRUE);
  111. }
  112.